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PREFACE 

This  thesis  defines  the  concept  of  a  numerical  system  for  a 
minicomputer  and  provides  a  description  of  the  software  and  computer 
system  configuration  necessary  to  implement  such  a  system.   A  procedure 
for  creating  a  numerical  system  from  a  FORTRAN  program  is  developed 
and  an  example  is  presented.   The  reader  should  have  some  knowledge  of 
FORTRAN  and  minicomputer  operating  systems .   A  familiarity  with  PAL 
assembly  language  for  the  PDP-11  is  necessary  to  make  full  use  of 
the  examples . 
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1.      THE  PROBLEM 

1.1  Definitions 

A  numerical  system  as  used  here  is  any  long  program  involving 
large  amounts  of  arithmetic  computation  on  a  set  of  numerical  data  that  is 
small  enough  that  manipulation  of  the  data  is  not  a  significant  part  of 
the  problem.   Manipulation  means  movement  of  the  data  from  one  storage 
device  to  another  or  change  from  one  format  to  another,  and  similar 
programming  problems  not  ordinarily  associated  with  numerical  analysis. 
Some  numerical  systems  would  be,  for  example,  programs  designed  to  solve 
difference  equations,  find  roots  of  nonlinear  equations,  perform  numerical 
quadrature  of  an  input  function,  minimize  a  function  of  several  variables, 
or,  in  the  example  dealt  with  here,  to  solve  a  system  of  ordinary 
differential  equations.   All  of  these  programs  have  in  common  the  require- 
ment of  a  great  deal  of  arithmetic  computation  and  simplicity  of  input: 
a  number  of  polynomial  coefficients,  a  difference  equation  of  specified 
order  and  coefficients,  a  subroutine  to  evaluate  a  function  with  some 
initial  conditions. 

1.2  Minicomputer  Use 

Such  a  numerical  system  is  often  difficult  to  write  for  a 
minicomputer  for  several  reasons.   Most  minicomputers  have  very  little 
computing  power  for  the  vast  amount  of  arithmetic  operations  needed  to 
run  a  program  of  this  type;  some  do  not  even  have  a  hardware  multiply/ 
divide  device,  making  complicated  arithmetic  operations  time  consuming. 
Small  word  sizes  of  12  or  16  bits  require  that  any  floating  point 
calculations  use  multi-word  manipulations  that  take  much  time  and  storage 


space.   The  small  size  of  main  memory  means  that  the  program  will  probably 
not  all  fit  into  it  at  one  time.  This  also  limits  the  amount  of  numerical 
data  since  storing  it  anywhere  hut  main  memory  would  be  impossibly  time 
consuming,  given  the  block  structure  and  slow  access  time  of  most 
peripheral  storage  devices. 

In  spite  of  these  drawbacks,  there  are  still  several  reasons 
for  creating  a  numerical  system  for  a  minicomputer.  Although  running  a 
numerical  system  is  time  consuming,  it  is  not  as  expensive  as  the  same 
amount  of  time  on  a  System/360  computer.   Some  minicomputers  are  idle 
part  of  the  day  and  a  set  of  numerical  systems,  possibly  in  a  program 
library  for  the  installation,  would  fill  in  the  gaps  by  being  run  in 
free  time  or  at  scheduled  times  daily.   Making  such  programs  available 
could  reduce  an  organization's  expenses  if  these  programs  are  normally 
taken  to  a  large  service  computer.   The  complexity  of  multi-word 
manipulation,  transcendental  function  evaluation,  and  programming  of 
complicated  formulae  is  reduced  greatly  because  many  l6-bit  word  size 
minicomputer  software  packages  include  both  a  library  of  subroutines  that 
perform  multi-word  manipulations,  and  a  FORTRAN  IV  compiler  that  allows 
execution  of  FORTRAN  programs  with  approximately  the  same  accuracy  as  a 
System/360  computer  using  single  or  double  precision  arithmetic. 

Finally,  the  problem  of  program  size  can  be  solved  by  proper  use 
of  peripheral  storage  devices  such  as  a  fixed  head  disc,  drum,  or  other 
"random"  access  device.  By  storing  large  program  segments  that  are  always 
executed  together  on  such  a  device  (hereafter  called  disc  although  any 
"random"  access  device 'may  be  substituted),  the  program  can  be  executed  by 
bringing  these  segments  into  a  core  buffer,  overlaying  the  program  segment 
previously  there.   This  core  buffer  will  be  called  the  overlay  buffer. 


1.3  Outline 

This  paper  will  deal  primarily  "with  a  procedure  for  transforming 
a  large  FORTRAN  numerical  program  known  to  work  properly  in  a  large  computer 
environment  into  a  numerical  system  that  runs  efficiently  in  a  minicomputer 
environment.   As  a  helpful  example,  a  problem  which  has  been  solved  using 
this  procedure  will  be  discussed  wherever  applicable.   The  coding  of  the 
example  numerical  system  is  found  in  Appendix  B  while  the  original  problem 
is  due  to  Gear  [3]  and  involves  integrating  a  first  order  system  of 
ordinary  differential  equations  through  some  interval  of  the  independent 
variable,  given  a  set  of  initial  conditions  and  a  subroutine  to  accomplish 
the  differentiation.   The  original  program,  DIFSUB,  is  designed  to  handle 
stiff  systems  of  equations  as  well  as  better  conditioned  ones. 

The  topics  to  be  presented  are  listed  below.   First,  a  computer 
system  which  is  well  adapted  to  the  problem  is  described;  criteria  for 
choosing  similar  systems  are  listed  as  well  as  references  to  several  such 
systems.   Next,  the  concept  of  a  program  overlay  is  developed.   The 
application  of  this  concept  in  partitioning  a  FORTRAN  program  will  be 
presented,  and  then  several  problems  resulting  from  such  a  partition  will 
be  discussed.  These  last  will  fully  describe  the  procedure  and  some 
additional  comments  about  it  will  conclude  the  paper. 


2.   THE  COMPUTER  SYSTEM  ; 

2.1  IGCS  Hardware  and  Software 

The  Illinois  Graphics  Computing  System  used  in  the  development 
of  DIF11,  the  numerical  system  described  in  Section  1.3,  belongs  to  the 
Department  of  General  Engineering  at  the  University  of  Illinois  at 
Urbana-Champaign .   As  described  in  [5]  and  [6],  it  is  based  on  a  PDP-11/20 
minicomputer  from  Digital  Equipment  Corporation  (DEC).   It  has  a  byte 
addressable  core  memory  of  l6K  l6-bit  words  with  a  ROM  bootstrap  loader. 
All  communication  between  the  CPU  and  peripheral  devices  takes  place  over 
a  single  bus  which  DEC  calls  UND3US.   Access  time  for  main  memory  is  490  ns. 
out  and  8  ns .  in.   Instructions  are  basically  of  the  two-address  type,  but 
12  different  addressing  modes,  eight  general -purpose  registers,  and  a 
software  stack  simulator  allow  a  wide  variety  of  instructions.   Instruction 
formats  require  one  to  three  words,  depending  on  the  addressing  mode  of 
operands.   A  KE11-A  extended  arithmetic  unit  provides  single  word  hardware 
multiply  and  divide  operations. 

A  256K  word  fixed  head  disc,  accessible  through  the  available 
software  in  64  word  blocks,  has  an  average  access  time  of  17  ms .  and  a 
transfer  rate  of  16  ys/word.   A  card  reader,  teletype,  and  Gould  4800 
electrostatic  printer /plotter  complete  the  basic  system.   The  Gould  4800 
is  used  by  DIF11  as  a  printer  through  a  software  program,  LP,  since  the 
Gould  will  eventually  be  the  output  device  of  a  system  containing  both  a 
CalComp-type  graphical  language  and  a  simulation  system  of  which  DIF11  will 
be  part  [5 ] • 

The  software  package  available  with  a  minicomputer  is  very  important 
since  it  can  greatly  reduce  the  work  needed  in  programming  a  large  system. 


IGCS  currently  operates  with  a  Disc  Operating  System,  DOS  version  UA, 
which  includes  utilities  for  transferring  data  "between  memory  and 
peripheral  devices,  an  assembler,  a  FORTRAN  IV  compiler,  and  a  relocating 
Linkage  Editor  which  presently  includes  an  overlay  facility.   See  [2] 
for  details.   EAELIB  is  a  library  of  FORTRAN  functions  (e.g.  SIN,  COS) 
and  utility  routines  (e.g.  floating  point  addition)  designed  to  operate 
with  the  KE11-A  multiply /divide  unit.   These  routines  are  linked  with  a 
FORTRAN  program  by  the  Linkage  Editor  and  are  used  by  the  program  at 
execution  time. 

This  system  is  suitable  for  a  numerical  system  since  l6K  of 
core  is  adequate  for  keeping  the  operating  system  core  (DOS  is  not 
entirely  resident,  but  overlays  part  of  itself  depending  on  its  current 
operation),  a  set  of  routines  from  EAELIB  for  a  large  program,  a  main 
program,  the  overlay  buffer(s),  and  storage  for  data  whose  size  is  variable 
and  dependent  on  an  input  parameter.   The  FORTRAN  compiler  is  very 
efficient  in  compiling  code  to  do  the  actual  computation  involved,  although 
control  transfer,  allocation  of  memory  space,  etc.,  are  best  performed 
using  assembly  language  interfaces.   FORTRAN  programs  interpret  a  polish 
string  of  addresses  of  routines  placed  in  core  by  the  Linkage  Editor  from 
EAELIB;  these  strings  are  the  result  of  compilation  of  the  program. 
Storage  is  best  handled  as  FORTRAN  arrays  whose  elements  fill  consecutive 
four  byte  segments  of  core  and  are  easily  processed  by  either  FORTRAN  or 
the  assembly  language  interfaces.   The  original  working  version  of  DIF11 
was  implemented  in  12K  core  with  user-written  overlay  software  so  the 
present  system  is  more  than  adequate.   The  speed  of  access,  relatively 
small  blocksize,  and  adequate  device  driver  software  make  the  fixed  head 
disc  easily  accessible  for  transferring  unchanging  program  segments  to  core, 


2.2  Alternate  Systems  - 

Inspection  of  the  above  system  leads  to  a  set  of  criteria  which 
a  minicomputer  system  should  satisfy  to  create  a  well  suited  environment 
for  a  numerical  system.   These  criteria  follow. 

The  main  memory  should  he  an  adequate  size  to  contain  the  system's 
resident  monitor,  a  large  section  of  code  from  the  numerical  system,  and  all 
the  variable  storage  for  the  system.   An  adequate  word  size  to  minimize  multi- 
word floating  point  manipulation  is  important;  at  least  16  bits  is  recommended. 
The  execution  speed  should  be  fast  enough  to  allow  an  adquate  turnaround  time 
for  the  size  problem  being  run;  a  good  criterion  would  be  that  the  minicomputer 
turnaround  time  including  loading  the  numerical  system  and  execution  should 
be  no  more  than  10  times  the  turnaround  for  a  full  size  service  computer  from 
input  to  retrieval  of  output. 

The  peripheral  device  that  acts  as  a  bulk  memory,  usually  a 
fixed  head  disc  or  drum,  should  have  fast  transfer  into  core.  Any 
sequential  device  would  be  unsuitable  for  this  . 

An  efficient  FORTRAN  compiler  (i.e.,  programs  compiled  use  close 
to  the  absolute  minimum  of  core)  is  necessary.   Since  many  numerical 
systems  may  need  to  use  double  precision  arithmetic,  this  feature  should 
compile  and  execute  efficiently. 

Finally,  a  good  software  monitor  capable  of  fast  and  effective 
control  of  device  transfers  while  occupying  a  minimum  amount  of  main 
memory  is  important,  both  at  execution  time  and  during  creation  of  the 
numerical  system. 

Several  computer  systems,  while  dissimilar  to  the  PDP-11/20 
in  many  ways,  appear  to  meet  the  above  criteria.   Some  may  be  better 


suited  to  the  purposes   of  the  purchaser  than   others   and  still  provide   an 
acceptable   environment   for  numerical   systems.      They  are: 

Hewlett-Packard  2100  A    [k] 

Interdata  Model  ?0    [7] 

Systems   72    [8] 

Varian  620/f    [9] 

Westinghouse  2500    [10] 


3.   THE  OVERLAY  CONCEPT 

3.1  Definitions 

A  program  overlay  is  a  segment  of  executable  code  intended  to 
be  used  with  a  permanently  resident  program  segment.   However,  this  code 
is  stored  permanently  only  on  a  peripheral  storage  device  (disc)  and  is 
brought  into  a  core  buffer  only  when  needed.   When  another  overlay  is 
needed  instead,  that  overlay  is  brought  from  disc  and  placed  into  the 
same  buffer,  thus  "overlaying"  the  previous  segment.   Clearly,  N  overlays 
of  approximately  equal  size  effectively  expand  the  overlay  buffer  to 
simulate  a  memory  area  N  times  as  large.  However,  one  pays  for  this  in 
terms  of  the  cost  of  disc  to  core  transfer  time.   This  is  more  advantageous 
than  having  all  executable  code  in  core  (if  possible)  and  using  the  disc 
to  store  variables  because  these  must  be  transferred  not  only  into  core, 
but  also  back  from  core  to  disc  after  being  processed,  whereas  overlays 
are  constantly  "refreshed"  from  an  unchanging  disc  image. 

Since  overlays  for  numerical  systems  are  usually  parts  of  a 
FORTRAN  program  originally  intended  to  be  executed  in  the  same  main  memory, 
the  overlay  should  have  access  to  all  variables  stored  in  resident  core, 
and  to  all  resident  subroutines,  especially  such  utility  routines  as  are 
normally  used  by  minicomputer  FORTRAN  compilers  for  complicated  arithmetic 
operations  such  as  floating  point  addition  (e.g.  EAELIB  in  Section  2.1). 
Because  of  these  features ,  certain  control  paths  between  the  resident 
code  and  each  overlay  must  be  provided.   These  paths  are: 

a.   Resident  to  overlay 

b  .   Overlay  to  overlay 

c.  Overlay  to  resident;  return  of  control 

d.  Overlay  using  resident  subroutines 


A  core  map  of  an  idealized  numerical  system  is  in 
Figure  3.1.1. 
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Figure  3.1.1.  Sample  Program  Map  with 
Control  Paths  Indicated 


3.2  User-written  Overlay  Facilities 

Some  minicomputer  software  systems  do  not  have  overlay  facilities; 
in  these  cases  the  user  must  code  his  own.   Such  a  system  was  written  for 
use  with  IGCS  before  versions  of  FORTRAN  and  the  Linkage  Editor  were 
provided  with  one.   The  code  for  this  system  in  PAL  assembly  language  is 
in  Appendix  A,  and  a  core  image  for  creation  of  an  overlay  file  with  it  is 
in  Figure  3.2.1. 
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Figure  3.2.1.   Core  Map  for  Creating  an 
Overlay  Using  MAKOVL 
The  system  consists  of  two  programs  called  CONCTM  and  MAKOVL. 
CONCTM  is  a  core  resident  subroutine  which  contains  tables  of  the  file 
names  of  all  overlays,  their  length  in  words,  and  their  starting  addresses 
in  core.   On  being  called,  it  looks  up  the  disc  start  block  of  the  overlay 
using  the  FILE  BLOCK  and  LINK  BLOCK  [2]  and  places  these  in  the  TRAN  BLOCK. 
The  overlays  are  contiguous  on  disc,  i.e.,  if  they  occupy  more  than  one 
disc  block,  they  are  placed  in  consecutive  disc  blocks.   When  an  overlay 
is  called,  the  call  is  to  CONCTM  with  the  first  argument  the  integer 
designating  the  overlay  wanted  and  the  following  arguments  being  the  proper 
arguments  for  the  subroutine  from  which  the  overlay  was  made.   CONCTM 
transfers  the  overlay  into  core  at  the  indicated  starting  address,  using 
a  TRAN  request  and  the  TRAN  BLOCK  mentioned  above.   Since  the  addresses  of 
the  arguments  of  the  overlay  subroutine  follow  the  statement  JSR  CONCTM, R5 
in  the  main  program,  and  since  subroutine  arguments  are  referenced  relative 
to  the  return  address  stored  in  register  R5,  a  simple  transfer  to  the  start 
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address  of  the  overlay  starts  its  execution.  A  FORTRAN  RETURN 
statement  shifts  control  back  to  the  resident  program  segment. 

The  overlays  are  created  "by  linking  together  all  resident 
parts  of  the  program,  the  overlay,  and  MAKOVL  below  these  (the  Linkage 
Editor  links  in  the  order  of  its  input  string,  filling  the  top  of  core). 
MAKOVL  is  executed  from  the  keyboard  and  computes  the  length  in  words 
and  disc  blocks  of  the  overlay  and  transfers  the  overlay  to  disc  under 
the  name  placed  in  its  FILE  BLOCK  (this  can  be  set  from  the  keyboard) . 
The  user  saves  the  length  and  starting  address  and  inserts  them  into  the 
final  version  of  CONCTM  when  the  routines  indicated  as  the  final  link  in 
Figure  3.2.1  are  linked  together  as  the  initial  load  module. 

Here  it  is  difficult  to  separate  the  floating  point  package 
routines  between  those  needed  by  the  resident  segment  and  those  needed 
only  by  the  overlays .   One  solution  is  to  determine  all  utility  routines 
needed  by  any  program  and  force  them  to  be  permanently  resident.   On  the 
PDP -11/20,  this  can  be  done  by  linking  each  program  and  inspecting  its 
load  map  which  lists  all  utility  routines  called  by  it.   The  union  of  all 
utility  routines  listed  on  any  load  map  is  then  placed  in  a  GLOBL  state- 
ment.  A  version  of  DIF11  using  the  above  system  was  programmed  and  tested. 

3.3  Manufacturer-supplied  Overlay  Facility 

The  overlay  facility  described  here  is  found  in  more  detail  in  [l], 
It  is  typical  of  those  available  from  manufacturers'  software  groups.   The 
system  uses  a  single  multiply  entrant  subroutine  LINK  with  entry  points 
LINK  and  RETURN.   These  provide  all  the  control  paths  described  in 
Section  3.1.   CALL  LINK('FILE')  in  resident  FORTRAN  code  or  JSR  LINK,R5  with 


12 

the  proper  argument  in  assembly  language  l)  initializes  FORTRAN  traceback 
routines,  2)  saves  register  0-5,  3)  causes  the  named  file  to  be  found, 
transferred  into  core,  and  executed  as  if  it  were  a  main  program.   The 
statement  CALL  RETURN  restores  register  0-5  and  returns  control  to  the 
place  in  resident  code  where  LINK  was  last  called. 

Since  it  is  easiest  to  write  and  test  the  FORTRAN  part  of  the 
numerical  system  with  the  overlays  as  subroutines,  a  problem  arises  in 
that  the  overlays  must  be  set  up  as  main  program  segments .   In  some  cases 
the  FORTRAN  SUBROUTINE  statement  can  be  removed,  but  if  any  parameters 
must  be  passed  as  if  the  overlay  were  a  subroutine  (see  Section  5«2), 
then  an  assembly  language  main  program  can  be  written  with  the  appropriate 
arguments.   Such  an  example  is  the  routine  INTRO  in  Appendix  B. 
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Figure   3.3.1.      Core  Map   for  Creating  an  Overlay 

Using  Manufacturer-supplied  Software 


The  overlays   are   created  using  the   relocating  Linkage  Editor, 
LINK.      The  resident  program  segments    are   first  linked  together,  bringing 


13 

all  floating  point  utilities  used  by  the  resident  code  into  core.   Two 
of  the  outputs  of  LINK  are  saved:   the  object  module,  and  the  Symbol 
Table  -which  lists  all  subroutine  start  addresses  including  those  for  the 
utility  routines.   Then  each  overlay  (with  its  main  program  INTRO  if 
needed)  is  created;  the  first  input  to  LINK  is  the  Symbol  Table,  and  the 
top  address  of  the  output  load  module  is  set  to  two  less  than  the  bottom 
of  the  core  resident  module  just  created.   The  Symbol  Table  input  allows 
the  overlay  to  use  utility  routines  already  core  resident  and  only 
utility  routines  which  it  alone  needs  are  added  to  the  overlay.   Figure  3.3.1 
shows  the  configurations  of  these  different  uses  of  LINK  to  create  the 
resident  and  overlay  modules . 
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h.      PARTITIONING  THE  FORTRAN  PROGRAM 

k.l     Criteria  for  Partitioning 

After  providing  for  an  overlay  facility,  one  next  considers  the 
overlays  themselves.   It  is  best  to  begin  -with  a  FORTRAN  program  which  is 
known  to  run  in  a  large  computer  environment  and  insure  that  the  overlay 
interconnections  will  work  by  breaking  the  program  down  into  overlay 
subroutines  and  running  it  again  on  the  same  computer. 

There  are  really  no  hard  and  fast  rules  that  can  be  applied, 
but  since  we  are  dealing  with  relatively  slow  minicomputers  and  one  of 
the  slower  operations  is  disc-to-core  transfers,  one  goal  of  the 
partitioning  process  is  to  minimize  the  number  of  such  transfers . 

Since  core  space  is  at  a  premium,  most  of  the  original  program 
should  be  put  into  overlays,  leaving  minor  bookkeeping  arithmetic  and 
major  control  transfer  decisions  in  the  resident  program  segment;  by 
making  the  overlays  of  nearly  uniform  length,  the  overlay  buffer  (which 
is  as  long  as  the  longest  overlay)  is  also  made  as  small  as  possible. 
The  implementations  of  these  two  goals  are  clearly  at  odds  since  the  real 
minimum  number  of  overlay  transfers  from  disk  would  maximize  the  overlay 
buffer,  and  the  minimum  overlay  size  would  be  only  two  or  three  FORTRAN 
statements,  requiring  innumerable  disc  transfers,  usually  of  length  less 
than  one  disc  block.  Proper  implementation  calls  for  a  balance  between 
these  two  extremes. 

U.2   COMMON  Variables 

Because  the  contents  of  the  overlay  buffer  at  the  time  a  new 
overlay  is  transferred  in  will  be  overwritten,  any  FORTRAN  variables 
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(other  than  special  constants)  must  be  core  resident.   This  is  most 
easily  accomplished  with  COMMON  statements  which,  when  linked  into  load 
modules  by  most  Linkage  Editor  facilities,  leave  all  of  the  variables 
in  a  resident  core  block.  When  the  overlay  modules  are  linked,  the 
addresses  of  elements  in  the  COMMON  block  are  available  from  the  Symbol 
Table  and  are  linked  with  the  overlay.   To  form  the  COMMON  block,  each 
overlay  and  the  resident  program  is  inspected  and  all  of  its  fixed  length 
variables  are  listed.   A  variable  array  is  of  fixed  length  if  its 
dimension  does  not  depend  on  an  input  to  the  program.   Next,  any  variable 
appearing  on  any  two  lists  must  obviously  be  shared  by  two  different 
routines  and  is  placed  in  the  COMMON  block.  Any  data  that  remains  the 
same  in  an  overlay  and  is  not  shared  by  any  other  overlay  can  be  defined 
by  a  DATA  statement  and  need  not  be  in  the  COMMON  block.   It  may  be 
convenient  to  have  several  COMMON  blocks  since  some  variables  may  be 
used  exclusively  by  certain  overlays  and  no  others  .  Allocation  of  space 
for  arrays  with  variable  dimensions  will  be  discussed  in  Section  5.3, 
but  note  that  the  initial  address  of  each  such  variable  must  be  passed 
to  the  overlay  as  a  subroutine  parameter  (it  cannot  be  in  a  COMMON  block 
since  its  initial  address  is  not  available  at  Linkage  time),  and  the 
variable  dimension  must  be  passed  to  the  overlay  as  a  subroutine  parameter 
as  required  by  most  FORTRAN  compilers . 

h.3     Program  Logic 

Some  few  FORTRAN  programs  are  executed  sequentially,  i.e.,  starting 
at  the  first  statement  and  "dropping  through"  to  the  end  with  a  few  tight 
DO-loops .   Partitioning  these  is  trivial.   The  program  is  divided  into  equal 
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length  overlay  segments;  either  each  overlay  calls  the  next,  or,  if 
this  is  not  possible,  a  short  resident  program  calls  each  overlay  in 
sequence.  However,  most  programs  are  not  this  easy. 

All  or  part  of  the  more  usual  program  is  executed  iteratively 
either  a  fixed  number  of  times  or  until  certain  conditions  are  met. 
DIF11,  for  example,  predicts  next  values  of  the  differential  variables, 
corrects  these  values,  and,  under  certain  conditions,  changes  other 
control  variables  (order,  step  size,  etc.).   This  whole  cycle — predict, 
correct,  change — is  repeated  until  a  final  value  is  reached  for  the 
independent  variable. 

The  best  approach  to  partitioning  such  a  program  is  to  follow 
a  principle  of  least  use;  if  something  is  seldom  used,  it  is  an  ideal 
candidate  for  an  overlay.   Many  programs  have  a  long  initialization 
section  which  is  only  executed  the  first  time  the  program  is  used. 
Clearly,  this  should  be  an  overlay.   The  original  version  of  DIF11  had 
such  an  overlay  but  the  present  version  does  not  need  one  since  the 
initialization  was  shortened.  A  segment  that  is  executed  at  most  once 
per  iteration  is  the  logical  next  candidate.   In  DIF11,  SUB3  and  SUB5  are 
executed  only  once  per  iteration  and  SUB1,  SUB2,  and  SUB6  are  executed 
only  when  a  change  of  control  variables  is  being  considered  or  carried 
out.   (See  Appendix  B). 

DO-loops  are  best  placed  completely  within  an  overlay  to  avoid 
each  pass  through  the  loop  causing  one  or  more  disc  transfers  to  the 
Overlay  buffer.   Occasionally,  however,  a  DO-loop  will  contain  enough 
material  for  several  overlays  (remembering  that  uniform  overlay  size  is 
a  goal  of  the  partition),  and  the  DO-loop  then  contains  several  overlays. 
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DIFFUN,  SUBT,  and  SUBU  are  all  contained  in  one  loop  in  DIF11.   The  same 
loop  contains  code  for  the  numerical  evaluation  of  the  Jacobian  matrix 
which  could  have  been  an  overlay.   It  was  not  because  this  would  have 
required  one  overlay  to  call  another  and  to  return  to  the  original 
overlay.   This  complicated  process  is  best  avoided  but  could  be  handled 
by  using  two  overlay  buffers.   See  Section  5.2. 

Occasionally,  one  especially  long  overlay  may  be  needed  but 
only  part  of  it  need  be  in  core  the  whole  time  it  is  being  executed.   In 
such  a  case  only  enough  of  the  overlay  should  be  brought  into  core  at 
one  time  to  fill  the  regular  sized  overlay  buffer,  and  when  some  part  of 
it  that  will  not  be  executed  again  is  completed,  the  rest  is  overlayed 
there.   This  will  be  possible  if  the  overlay  facility  allows  overlays  to 
call  overlays  as  specified  in  Section  3.1. 

Any  segment  that  is  called  from  a  number  of  different  places  in 
the  resident  code  is  probably  best  left  in  core,  but  if  it  is  long  or 
called  very  seldom,  it  can  be  made  an  overlay  as  in  SUB8  in  DIF11. 

k.h     Control  Transfer 

It  must  be  noted  that  each  overlay  may  have  several  GO  TO 
statements  which  reference  a  statement  in  the  resident  module,  or  worse, 
in  another  overlay.   The  partition  should  minimize  this,  but  if  it  happens, 
than  a  control  transfer  vector  can  be  used  to  allow  this.   Either  a  single 
variable,  JLINK,  or  an  array,  JLNK(N),  can  be  used;  the  various  values 
should  have  some  significance  to  the  program.   The  single  variable  is 
probably  better  since  it  is  less  complicated.   Set  JLNK  =  1  on  entering 
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an  overlay,  then  if  a  statement  is,  say 

GO  TO  750 
and  750  is  in  the  resident  module,  change  the  statement  to 

JLNK  =  2 
RETURN 
If  a  statement  is 

IF( LOGICAL  EXPRESSION) GO  TO  715 
change  this  to 

IF (LOGICAL  EXPRESSION )JLNK= 3 
IF ( JLNK. EQ. 3) RETURN 
This  requires  that  any  arithmetic  IF  statements  he  made  logical  IF  statements 

Immediately  after  the  overlay  calling  statement  in  the  resident 
code,  a  computed  GO  TO  statement  can  he  used  to  effect  the  transfer: 

CALL  OVERLAY  (SUB 5) 
GO  TO(5140,750,715,..  .),JLNK 
5^0   ... 
remembering  the  possibility  of  executing  the  next  sequential  statement  if 
JLNK  remains  1. 

Transfers  into  the  middle  of  another  overlay  can  be  handled  by 
a  computed  GO  TO  at  the  beginning  of  the  overlay  so  that  this  is  the 
first  statement  other  than  subroutine  calls  after  the  computed  GO  TO  in 
the  resident  code.   JLNK  must  be  set  to  1  after  executing  the  initial 
computed  GO  TO  statement.   For  example, 


GO  TO  TOO 
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"becomes 


JLNK=2 
RETURN 
*  *  *  _  ********** 

5U0   GO  TO(550,600),JLNK 

600   CALL  OVERLAY (SUB2) 
***   ********** 


( overlay  1 ) 


(MAIN) 


GO  TO(650,TOO) ,JLNK 


TOO  JLNK=1 


(overlay  2) 


Examples  of  all  these  techniques  are  in  SUB5  of  DIF11  in  Appendix  B 
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5.   USER  SUPPLIED  OVERLAYS 

5.1  Us er  Input 

Many  numerical  systems  can  operate  with  a  minimum  of  input. 
For  example,  the  only  input  needed  to  solve  a  difference  equation  is  the 
order,  the  coefficients,  and  the  initial  values.   These  can  all  be  input 
on  cards  or  keyboard,  and  the  storage  for  this  data  can  be  in  the  core 
resident  main  program  by  specifying  a  maximum  size  for  each  variable 
array.   Other  programs  will  require  more  from  the  user  than  numerical 
data;  a  numerical  quadrature  program  will  require  that  some  way  be 
provided  to  evaluate  the  function  at  any  given  point,  usually  a  subroutine 
Some  programs  will  have  large  data  arrays  which  differ  greatly  in  size 
depending  on  input  parameters .   Storing  these  within  the  main  program 
segment  may  be  too  restrictive.   An  LU  matrix  decomposition  would  need 
such  a  variable  sized  space  to  store  the  matrix  and  outputs.   And,  some 
systems  require  both  an  evaluating  subroutine  and  variable  sized 
arrays,  as  is  the  case  with  DIF11. 

5.2  User  Subroutines 


The  solutions  to  these  problems  for  DIF11  are  presented  here; 
they  seem  applicable  to  similar  cases  on  other  systems  and  are  at  least 
a  good  starting  point.   In  DIF11,  the  user  is  responsible  for  writing  a 
FORTPAN  subroutine  called  DIFFUN(T, Y,DY) ,  and  by  following  carefully 
detailed  instructions,  creating  an  overlay  from  it.   First,  compile 
DIFFUW(T,Y,DY)  where  T  is  the  independent  variable,  DY(l)  is  the  first 
derivative  of  Y(l),  and  all  derivatives  of  order  m  higher  than  one  have 
been  replaced  by  m  first-order  equations  using  standard  techniques. 
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Y(l)  and  DY(l)  should  be  dimensioned  N,  the  number  of  equations.   Next,  call 
LINK,  the  Linkage  Editor.   Specify  DIFFUN. LDA  as  output;  specify  as  input 
ST  (the  Symbol  Table  of  Section  3.3  which  is  kept  available  on  disc  or 
tape),  the  object  module  of  INTRO  (also  available),  the  object  module 
resulting  from  compiling  DIFFUN,  and  the  EAELIB.  Also,  specify  that  the 
top  core  address  available  to  this  module  is  two  less  than  the  bottom  of 
the  resident  code  (this  address  is  posted  in  the  machine  room).   After 
running  LINK,  run  the  program  DIF11  with  two  data  cards:   one  specifies 
N,  T  ,  H  (the  length  of  the  integration),  EPS  (the  error  criterion), 
and  ISTORE  (the  octal  top  address  to  be  filled  by  dynamically  allocated 
variable  arrays);  the  other  lists  the  N  initial  values  of  Y(l). 

This  overlay  uses  the  same  overlay  buffer  as  all  other  overlays . 
This  is  possible  because  none  of  the  other  overlays  call  it.   Any  other 
overlay  could  call  DIFFUN  and  then  return  to  the  main  program;  however, 
a  great  deal  more  work  would  be  required  to  let  an  overlay  call  DIFFUN, 
have  DIFFUN  reinstate  that  overlay,  and  then  transfer  to  the  statement 
after  the  call  to  DIFFUN.   This  would  be  equivalent  to  a  regular  subroutine 
call  to  DIFFUN. 

Another  possibility  is  to  provide  a  second  buffer  just  for  the 
user  supplied  overlay,  allowing  all  overlays  to  call  it.   However,  it 
must  be  remembered  that  this  takes  valuable  core  that  might  otherwise  be 
used  for  variable  storage.   Since  there  is  no  way  to  predict  the  length 
of  a  user  supplied  module,  this  allows  the  first  overlay  buffer  to  be  of 
fixed  length  while  the  second  varies.   This  may  be  an  important  consideration 
for  some  users. 
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5.3  Variable  Size  Arrays 

The  problem  of  variables  whose  dimensions  depend  on  an  input 
parameter  is  best  solved  by  allooating  storage  space  for  the  FORTRAN 
arrays  in  an  assembly  language  MAIN  program.   In  DIF11,  various  arrays 
are  dimensioned  (N,T),  (N,2),  or  (N,N) .   The  MAIN  program  (see 
Appendix  B)  computes  U*N  (each  floating  point  -word  uses  four  bytes  on 
the  PDP-11/20),  28*N,  8*N,  and  U*N**2.   Using  RO  (register  0)  initially 
set  to  ISTORE  as  a  pointer,  MAIN  subtracts  an  appropriate  number  of 
bytes  from  RO  and  assigns  RO  to  an  address  vector  ADDR .   Then  it  places 
the  addresses  of  the  appropriate  variables  in  the  address  list  of  the 
calling  sequences  of  INPUT,  OUTPUT,  and  DIFSUB— the  FORTRAN  subroutines 
called  from  MAIN — and  executes.   It  should  be  noted  that  most  FORTRAN 
compilers  require  a  variably  dimensioned  array,  e.g.  Y(N,7),  to  be  in 
a  subroutine  with  the  variable  dimension  N  as  an  argument,  so  DIFSUB, 
INPUT,  and  the  overlays  must  be  subroutines.   Whenever  DIFSUB  calls  one 
of  its  overlays ,  the  overlay  starts  with  a  main  program  INTRO  as 
required  by  the  overlay  facility,  and  INTRO  places  the  initial  addresses 
assigned  to  the  variably  dimensioned  arrays  and  kept  in  an  assembly 
language  named  .CSECT  block  (same  as  FORTRAN  COMMON)  into  the  FORTRAN 
calling  sequence.   Then  INTRO  executes  the  subroutine  call,  returning  to 
the  overlay  main  program  before  calling  RETURN  (see  Section  3.3)  in 
order  to  restore  all  registers.   In  this  way,  eight  variably  dimensioned 
arrays  are  allocated  dynamically  in  a  manner  similar  to  ALGOL  by  using 
assembly  language  interfaces  between  FORTRAN  program  segments . 

It  should  be  noted  that  the  pseudo-stack  mentioned  in  Section  2.1 
is  kept  immediately  below  the  overlay  buffer  by  the  DEC  provided  software, 
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and  extends  into  lower  core  addresses  .   This  requires  that  the  dynamically- 
allocated  variables  he  placed  well  below  the  stack  and  is  the  reason 
ISTORE  is  an  input  parameter.   If  a  user  written  overlay  facility  were 
available,  the  stack  could  be  moved  at  will  below  the  data. 
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6.   CONCLUDING  OBSERVATIONS 

The  following  observations  seem  appropriate  to  the  problem 
discussed  and  are  presented  in  an  attempt  to  provide  some  perspective. 
Firstly,  the  method  of  reducing  a  large  FORTRAN  program  to  a  mini- 
computer overlay  system  is  best  considered  as  a  procedure.   A  list  of 
step  by  step  instructions  can  be  found  in  the  preceding  chapters ,  but 
no  guarantee  is  provided  that  the  procedure  will  ever  stop.   Thus, 
this  is  not  an  algorithm  with  input:   one  large  FORTRAN  program; 
output:   equivalent  minicomputer  numerical  system.   A  counter  example 
is  a  program  with  larger  storage  requirements  than  could  ever  fit  into 
core  or  be  profitably  allocated  to  disc.   This  procedure  can  be  helpful 
but  no  promises  are  made . 

Secondly,  in  most  cases  the  more  general  a  program  is,  the 
larger  it  becomes.  A  simple  Runge-Kutta  (R-K)  program  with  step 
doubling  will  take  much  longer  to  integrate  a  complicated  system  on  a 
360/75  than  DIFSUB  will;  but  if  a  problem  is  fairly  simple,  a  R-K 
program  requiring  no  overlays  would  be  much  faster  than  DIF11  on  the  PDP -11/20, 
Thus,  one  can  save  time  by  having  available  a  simple  method  for  simple 
problems  as  well  as  a  complicated  method  to  solve  complicated  problems. 

Finally,  one  should  be  wary  of  manufacturer  supplied  overlay 
facilities,  especially  if  overlays  are  called  more  than  a  few  times.   The 
overlay  processor  supplied  by  DEC  looks  up  the  current  disc  start  block 
of  each  overlay  every  time  it  is  called,  even  though  some  overlays  are 
called  hundreds  of  times'  from  the  same  disc  location.   For  this  reason, 
the  final  simulation  language  package,  ILLISM,  which  DIF11  will  be  part 
of,  will  use  a  locally  written  overlay  facility  that  looks  up  disc  addresses 
once,  stores  them  in  a  directory,  and  thus  executes  twice  as  fast. 
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With  these  comments  in  mind,  one  can  conclude  that  many  large 
numerical  analysis  programs  can  he  profitably  used  on  a  minicomputer 
system  using  the  techniques  described  here. 
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APPENDIX  A 
User-written  Overlay  Facility 
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QNO  NUMBER  OF  WORDS. 


PSR   NUMSEQ 
fiSR    NUMSEQ 
ODD   #1 , NUMSEQ 
ALLOC   LNKBLK, FILBLK.NUMSEG 
MOV    NUMSEQ, -<SP> 
MOV    #FILBLK,-<SP> 

MOV    #LNKBLK, ~<SP> 

EMT  1  5 

INC   <SP>  + 

BN  E   ERR 

.OPENC   LNKBLK,  F1L9I.K 


JNE   ERR 

.OPENC   LNKBLK, Fll 
MOV   #13,FILBLK-2 
MOV  #FILBLK, -<SP> 
MOV   #LNKBLK, -<SP> 
EMT  16 

MOV   #BUF,-<SP> 
MOV  #LNKBLK,-<SP> 
EMT  2 
WRIT   LNKBLK 
MOV   #LNKBLK, -<SP> 
E"  M  T  1 

.CLOSE    LNKBLK 
MOV   *LNKBLK, -(8P> 
EMT  1  7 


ERR: 

OVRLEN:     .WOKD 
NUMSEQ:    .WORD  0 
L  N  N K  BLOCK 
•WORD  0 


.RLSE   LNKBLK 
MOV   #LNKBLK, -<SP> 
EMT  7 
.EXIT 
EMT  60 
.  WORD 
0 


.WORD  0,0,. 
•  DF/ 


1 


LNKBLK: 

.RPD50   ' vr  r 
;  FILE  BLOCK 

. WORD  0, 0 
F1I.BLK:    .RRD5B   'DIF' 

.RRD50   •FUN/' 
.RAD50  'LDA/ 

.BVTE   0,0,0,0 
BUF:       .WORD   OVRLEN, 7, OVRLEN, I 

.  END 


'IFFUN 
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SUBROUTINE  CONCTM 
•TITLE   CONCTM 

SP  «  16 


CONCTM: 


POINTR: 


SP  «  *6 

R0  "   f e 

GLOBL  CONCTM 

S2<*5> .POINTR 
POINTR 
POINTR 
R0, -<SP> 


MOV 
OSL 
ASL 
MOV 
MOV 
MOV 
MOV 
MOV 
MOV 
MOV 
MOV 
EMT 
MOV 
CLR 
MOV 
EMT 
MOV 
TST 
TST 
MOV 
MOV 
EMT 
MOV 
EMT 
MOV 
EMT 
JMP 


POINTR, R8  5  PUT  9EC0I 
POINTR<R0> , FILBLK+2  JFJLE  Bl 
P0INTR-2<R8>, TILBLK  ;  PUT 
ADDR<R8>  ,  TRNBLK*4 
ADDR-2CR0)  ,  TRNBLK+2  {PUT 
<SP>+,R0 


JMOVE  Sft 

jPUT  4*POINTR 

WORD  OF 


ON 


INDICATOR  TO  POINTR 
THE  STACK. 
MODULE  5045 
IN  FILE  BLOCK 

TRAN  BLOCK 


SECOND 
"  BLOCK 

FIR9T  WORD 


( SP)  ♦ , R0 
iLNKBLK, -<SP) 

A. 

-<SP) 
-<SP> 

, TRNBLK 


STORT  LOCATION  IN 


•FILBLK 

-<SP) 

•LNKBLK 
1  4 

<SP> 
<SP>* 
<SP)  + 


.  WORD 
, RAD50 
. RAD50 
.ROD50 
. RHD50 
. RQD30 
. RQD90 
.  RAD50 
.  R  h  C>  5  e 
.  RAD50 
.  R  ft  D  5  O 
.  R0D50 
. RAD50 
. RAD50 
. ROD50 
.  R  H  D  5  * 
R  ft  0  ">  0 
NUMWRD:ADDR: 

.  WORD 
.  WORD 
LNKBLK:     .WORD 
. RRD30 
.  WORD 
F ILBLK:     . WORD 
. RAD50 
.  B  V  T  E 
TRNBLK:     .WORD 
.  END 


#TRN8LK, -<SP> 
#LNKBLK, -<8P> 
1  0 

ttLNKBLK,  -<SP> 
1 

*LNKBLK, -<  SP) 
t 

STRNBLK+2 

0         ;PO 
y  D I F  / 

/FUN/ 

/PAS/ 

/COL/ 

/QFI/ 

/ND/ 

/SCfi/ 

/LE/ 

/EVfi/ 

/L/ 

/CHB/ 

/NGE  ' 

/INI/ 

/T/ 


.  I  N  IT 

;  PUT 

.WRITE 
.  WRIT 
.  ALSE 


•  LOOK 
START 


BLOCK  IN  TRAN  BLOCK 


5    .RLSE 

J    00  TO  SUBROUTINE  STORT  0 
INTER  BND  LIST  OF  FILE  NAMES 


ADDRESS 


/COR/ 

/ECT/ 
WORD  0, 37260, 76,37122, 155,371 12, 161 ,36476, 
36634,318,3  9064, 1 174,37072, 161 ,36414,420 

0      ;link  block 

0,0,  1 
/DF/ 
0,4      ; FILE  BLOCK 


367 


0,0.  1 
/DF/ 
0,4 
0,  0 

/LDA/ 
0,0,0,0 
0,0,0,4 
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APPENDIX  B 
DIF11  Program  Code 
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THIS  PROGRAM  IS  THE  MRIN  PROGRAM  OF  0 1 F t  1  . 
IT  INITIALIZES  CONSTANTS,  READS  IN  VARIABLES, 
AND  ASSIGNS  STORAGE  LOCATIONS  TO  VARIABLY 
DIMENSIONED  STORAGE  FOR  DIFSUB. 


.TITLE  MAIN 

Re  -  se 


R1 
R2 
R3 
R4 

R5 
SP 
PC 


f  1 
*2 

I  3 
*4 
*5 
X6 


.GLOBL 
. CSECT 
BEGIN: 


LOW), L0W2, OUT PUT, INPUT, LINK, RE TURN, N, DIFSUB 


.WORD  a 
RET1  : 


ON: 


NORD  0 


LISTO: 
RET2: 

LOOPS: 


L00P4 


ST  1  9: 
ST20: 


LIST: 
out: 

OUTLIS: 
RET3: 


N 
, INPUT 
RET1 
,  L.0W2,  XLIM 

LOW1 , L0W2 

ON 

L0W2, LOW  I 

N,  R1 

R1 

R1 

LOH1 , R0 

•  2,  RO 
R  I  ,  R0 
R8, ADV 
Rt  ,  R0 
RB.AERROR 

UB  R  1  , R0 
R0.AVMAX 
R1  ,  R0 
RB,  AI 
R1  ,  R3 
R3 
R3 
R3 

R1  ,  R3 
R3,  R0 
R8,  AV 
R3,  R0 
R0, ASQVE 
R1 

R1  ,  R0 
R0, ACSAVE 
N2,  R0 
R9, APM 

♦  4 


PIP 


CLR 
SR  R5 

BR 
,N,N2 

CMP 

BMI 

MOV 

MOV 

ASL 

ASL 

MOV 

SUB 

SUB 

MOV 

SUB 

MOV 
S 

MOV 

SUB 

MOV 

MOV 

ASL 

ASL 

ASL 

SUB 

SUB 

MOV 

SUB 

MOV 

ASL 

SUB 

MOV 

SUB 

MOV 

BR 

MOV  AV, LISTO 

JSR  R3, INPUT 

ER  RET2 

. WORD  0, N, 0, 0, 0 

MOV   AVMAX.R1 

MOV  N,R0 

MOV  #40200,  CR1  )♦ 

CLR  <R1 >♦ 


CLR  ( R1 

SUB  # 1 , R0 

CMP  #0,R0 

BNE  LOOPS 

MOV  AV, OUTLIS 

MOV  #20, R3 

MOV  #ADDR,R0 

MOV  #LIST,R1 

MOV  (RB)  +  ,  (R1  >■» 

TCT  CDOS- 


L  U  U  K 

MOV  AV.OUILI 

MOV  #20, R3 

10V  #ADDR,R0 

.10  V  •  L  I  S  T  " 

MOV  <R8>* 

TST  <R3>- 

C  M  P  #  0  ,  R  3 

BNE  L00P4 

CLR  JSTART 

BR  ST20 

MOV  #1 , JSTART 

JSR  R5, DIFSUB 

BR  OUT 

.WORD  N 

. WORD  0, 9, *, 8, 0, e, 

JSR  R3, OUTPUT 

BR  RET3 

. WORD  0, N, XLIM 

CMP  N  t  , KFLAG 

BEQ  3T1 9 

r  u  t  /.a 


;CALL  INPUT  SUBROUTINE  WITH  N-0  TO  READ  IN 
JN, T.H.EPS, AND  L0W2-BOTTOM  ADDRESS  OF 
{DIFFUN  OVERLAY  PROVIDED  BV  USER. 

[LOWEST  OVERLAY  START  ADDRESS  GOES  TO  L0W1 
JTO  BECOME  TOP  OF  STORAGE  AREA. 

{COMPUTE  ADDRESS  OF  DYNAMICALLY  ALLOCATED 
{STORAGE  FOR  DIMENSIONED  VARIABLES. 

;R1  CONTAINS  4*N, R0CONTAINS  TOP  OF  STORAGE. 

JADV  -  •  OF  DY<N> 

;AERROR  -  #  OF  ERROR<N> 
JAVMAX  «■  •  OF  VMAX(N> 
j  AIP  •  •  OF  IP<N> 


',  R3  CONTAINS  7  +  4*N 

; AV  -  #  OF  Y<N, 7> 

;ASRVE  •  •  OF  SAVE(N,7) 
;R1  CONTAINS  0*N 

{ACSAVE  •  #  OF  CSAVE<N,2) 

; N2  -  (N*+2>*4 

JAPW  -  «  OF  PW(N,N) 


;  READ  V(*,I)  FROM  CARDS 
;  SET  VMAX<*>  -  1.0 

JR0  18  THE  LOOP  COUNTER 

{PREPARE  PARAMETER  LISTS  FOR  CALLING 

{DIFSUB  AND  OUTPUT 


{ JSTART 
{ JSTART 


6,0 
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N-      .word  e,0 

N2:        • WORD  8, 8 

LOM1 :      . WORD  03351 2,0 

L0W2:      .WORD  0,0 

XLIM:      .WORD  0,0 

.CSECT  COM 

JSTPRT  -  .+40 

KFLOG  ■  .+50 

T  =  .+74 

H  =  .+124 

HMHX  -  .+140 

HMIN  »  .+144 

.CSECT  flDDR 

HC'DR:flDV:     .WORD  0 

PERROR:   .WORD  0 

flVMflX:     .WORD  0 

BIP:       . WORD  0 

a  v :  .  w  0  R  D  Q 
pspve:  .word  0 
pcspve:   .word  0 

PPW:  .WORD    0 

.END     BEGIN 


Resident  I/O  Routines 


.T.TOLD, 
OWN, ENQ1 ■ 


FlflQ.NQ.NQOLD.NEWQ, EPS, T, TOLD, 
HNEW,HMQX,HMIN,E,EUP,EDWN,ENG>1, 


Resident  FORTRAN  Main  Program 


33 


r 

SUBROUTINE  DIFSUB<N,DV, EAR OR, VMAX, IP, V, SAVE, C SAVE, PW> 

COMMON  /COM/-  A<7>,JLNK, JSTART, K,KFLAG,NQ,NQOLD,NEWC!,EPS,T, TOLD, 

1    IDOUB,  IWEVAL,  I  RET  ,  IRET1  ,H, HOLD, HNEW.HMAX, HMIN, E.EUP.EDWN, EN  C'1  , 

2ENQ2,ENQ3,BND,BR,DEL,DEL1  ,  D 

r+*»*  +  +  +  +  +  *  +  +  *  +  +  +  +  +  ***  +  **  +  +  ****  +  *  +  *  +  +  *  +  *  +  *  +  *4-***  +  +  +  +  **  +  *±**  +  +  +  **  +  +  .t.  +  +  4-** 

c  *  * 

C*  THIS  SUBROUTINE  INTEGRATES  A  SET  OF  N  ORDINARY  DIFFERENTIAL  FIRST  * 

C*  ORDER  EQUATIONS  OVER  ONE  STEP  OF  LENGTH  H  AT  EACH  CALL.    H  MAV  BE  + 
C*  INCREASED  OR  DECREASED  WITHIN  THE  RANGE  HMIN  TO  HMAX  TO 

C*  ACHIEVE  AS  LARGE  A  STEP  AS  POSSIBLE  WHILE  NOT  COMMITTING  A  SINGLE  + 

C  +  STEP  ERROR  WHICH  IS  LARGER  THAN  EPS  IN  THE  L-2  NORM,  WHERE  EACH  + 

C*  COMPONENT  OF  THE  ERROR  IS  DIVIDED  BV  THE  COMPONENTS  OF  VMAX.  + 


C  + 

c* 


THE  PROGRAM  REQUIRES  THE   SUBROUTINES  NAMED  * 

C*DIFFUN<T,V,DY>  * 

C*  DECOMP(N,M,PW,IP)  * 

C *  SOLVE'N  H  PW  CSAVEd  I)  IP)  * 

C*  THE  FIRST, 'dIFFUN,  EVALUATES  THE  DERIVATIVES  OF  THE  DEPENDENT  + 

C*  VARIABLES  STORED  IN  V  <  1  ,  I  >  FOR  I  -  1  TO  N,  AND  STORES  THE  * 

C+  DERIVATIVES  IN  THE  ARRAV  DV.  DECOMP  IS  A  * 

C*  STANDARD  LU  DECOMPOSITION  WITH  PIVOTING  THAT  DECOMPOSES  THE  MATRIX   * 

C*  PW,   LEAVING  THE  PIVOTS  IN  THE  INTEGER  ARRAV  IP.    M  IS  THE  DECLARED   * 

C*  SIZE  OF  PW.    IP<N>  IS  SET  TO  0  IF  PW  IS  SINGULAR.    SOLVE  PERFORMS  * 

C  *  BACK  SUBSTITUTION  ON  THE  CONTENTS  OF  C  S  A  V  E  <  I  ,  1  >  ,  LEAVING  THE  * 

C*  RESULTS  THERE.  + 

C  *  * 

C*  THE  TEMPORARV  STORAGE  SPACE  IS  PROVIDED  BV    *MAIN*     IN  THE  * 

C+  INTEGER  ARRAV  IP,   THE  SINGLE  PRECISION  ARRAVS  PW,   DV,  + 

C*  PRECISION  ARRAVS  SAVE  AND  C8AVE.    THE  ARRAV  PW  IS  USED  ONLV  TO  HOLD  + 

C*  THE  MATRIX  OF  THE  SAME  NAME,  AND  SAVE  IS  USED  TO  SAVE  THE  VALUES  * 

C*  OF  V  IN  CASE  A  STEP  HAS  TO  BE  REPEATED,  BUT  CSAVE  IS  USED  TO  HOLD  * 

C*  SEVERAL  ARRAVS.  * 

C  +  * 

C*  THE  PARAMETERS  TO  THE  SUBROUTINE  DIFSUB  HAVE  * 

O  THE  FOLLOWING  MEANINGS..  + 

C+  * 

C*  N         THE  NUMBER  OF  FIRST  ORDER  DIFFERENTIAL  EQUATIONS.    N  + 

C*  MAV  BE  DECREASED  ON  LATER  CALLS  IF  THE  NUMBER  OF  * 

C  *  ACTIVE  EQUATIONS  REDUCES,  BUT  IT  MUST  NOT  BE  * 

C*  INCREASED  WITHOUT  CALLING  WITH  JSTART  -0.  * 

C*  T         THEINDEPENDENTVARIABLE.  * 

C*  V          A   7  BV  N  ARRAV  CONTAINING  THE  DEPENDENT  VARIABLES  AND  * 

C *  THEIR  SCALED  DERIVATIVES.    Y  <  J  ♦  1  ,  I  >  CONTAINS  * 

C*  THE  J-TH  DERIVATIVE  OF  V<I>  SCALED  BV  * 

C  *  H  +  +  j'FACTGRIAL<'J>  WHERE  HIS  THE  CURRENT  * 

C  *  STEP  SIZE.  ONLV  V<1,I>  NEED  BE  PROVIDED  BV  * 

C  +  THE  CALLING  PROGRAM  ON  THE  FIRST  ENTRV.  * 

C  *  H          THE  STEP  SIZE  TO  BE  ATTEMPTED  ON  THE  NEXT  STEP.  + 

C*  H  MAV  BE  AD.  JUSTED  UP  OR  DOWN  BV  THE  PRO GRAM  * 

C*  IN  ORDER  TO  ACHEIVE  AN  ECONOMICAL  INTEGRATION.  * 

C*  HOWEVER,   IF  THE  H  PROVIDED  BV   ♦MAIN*   DOES  * 

C*  NOT  CAUSE  A  LARGER  ERROR  THAN  REQUESTED,   IT  * 

'"*  WILL  BE  USED.    TO  SAVE  COMPUTER  TIME,  *  M  A  I  N  *  * 

C*  USESPFAIRLVSMALLSTEPFORTHEFIRST  * 

C  *  CALL.    IT  WILL  BE  AUTOMATICALLY  INCREASED  LATER.  * 

C*  HMIN      THE  MINIMUM  STEP  SI2E  THAT  WILL  BE  USED.  * 

C*  UMAX      THE  MAXIMUM  SIZE  TO  WHICH  THE  STEP  WILL  BE  INCREASED  * 

C*  EPS       THE  ERROR  TEST  CONSTANT.    SINGLE  STEP  ERROR  ESTIMATES  * 

C  *  D  I  V  I  D  E  D  B  V   Y  M  A  X  <  I  >    MUSTBELESSTHANTHIS  * 

C  *  IN  THE  EUCLIDEAN  NORM.  THE  STEP  AND'OR  ORDER  IS  * 

C*  ADJUSTED  TO  ACHEIVE  THIS.  * 

C  *  VMAX      AN  ARRAV  OF   N  LOCATIONS  WHICH  CONTAINS  THE  MAXIMUM  * 

C*  OF  EACH  V  SEEN  SO  FAR.    IT  SHOULD  NORMALLY  BE  SET  TO  * 

C  *  1   IN  EACH  COMPONENT  BEFORE  THE  FIRST  ENTRV.   'SEETHE  * 

C*  DESCRIPTION  Or  EPb. ^  + 

C  *  ERROR     AN  ARRAV  OF   N  ELEMENTS  WHICH  CONTAINS  THE  ESTIMATED  * 

C*  ONE  STEP  ERROR  IN  EACH  COMPONENT.  * 

C  *  KFLAG    A  COMPLETION  CODE  WITH  THE  FOLLOWING  MEANINGS..  * 

C*  +1    THE  STEP  WAS  SUCCESFUL.  * 

C*  -1    UNRECOVERABLE  ERROR  + 

C*  0    INTEGRATION  COMPLETED  * 

C*  JSTART   AN  INPUT  INDICATOR  WITH  THE  FOLLOWING  MEANINGS..  + 

C*  0   PERFORM  THE  FIRST  STEP.    THE  FIRST  STEP  * 

C*  MUST  BE  DONE  WITH  THIS  VALUE  OF  JSTART  * 

C*  SO  THAT  THE  SUBROUTINE  CAN  INITIALIZE  * 

C*  ITSELF.  * 

C*  *1    TAKE  A  NEW  STEP  CONTINUING  FROM  THE  LAST.  * 

C*  JSTART  16  SET  TO  NQ,  THE  CURRENT  ORDER  OF  THE  METHOD  * 

C+  AT  EXIT.    NQ  IS  ALSO  THE  ORDER  OF  THE  MAXIMUM  * 

C*  DERIVATIVE  AVAILABLE.  * 

C  *  PW        A  BLOCK  OF  AT  LEAST  N**2  FLOATING  POINT  LOCATIONS.  * 

i  *********  +  **♦***  +  ****♦***  +  ♦  +  **♦****  +  *  +  +  ****  +  **  +  *  +  +  +  **♦  +  +  +  +  ****  +  +  +  .♦.**  +  ♦* 

DIMENSION  V<N,7>,VMAX(N>,SAVE<N,7),ERR0R<N),PW<1  >,DY<N># 
1    CSAVECN, 2) , IP(N) 
IRET  -  1 
KFLAO  ■  1 
IF  k JSTART. LE.0>  GO  TO  140 
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DO  110  I  »  1 ,N 
DO  lie  J  -  1  ,K 

SAVEC I , J)  -  V( I , J) 
HOLD  -  HNEW 

IF  '.  H  .  E  Q  .  H  0  L  D  )  GO  TO  130 
120  I  RET  1  -  1 
GO  TO  7  5  0 
I  3  6   N  0  OLD  =  N  Q 
TOLD  =  T 

IF  (JSTART. 6T. 8)  GO  TO  250 
GO  TO  170 
140   IF  (JSTART.EQ.-I)  GO  TO  160 
r*  +  *******  +  +  *  +  **  +  *  +  ***  +  *  +  ***  +  *  +  *  +  *****#  +  ********** +  +  ****  +  **  +  +  +  *  +  +  +  **  +  +  ■*•+ 
O  ON  THE  FIR3T  CALL,  THE  ORDER  IS  SET  TO  1  AND  THE  INITIAL  + 

C*  DERIVATIVES  PRE  CALCULATED.  * 

BR  -  1.0 
HO.     -  1 
C       CALL  DIFFUNCT, V, DV) 
CALL  LINK<  '  DIFFUN  '  > 
DO  150  I  -  1,N 
1 50    V< I , 2)  -  DV< I )*H 
HNEW  -  H 
K  -  2 
GO  TO  100 


BO   TO   100 

C*  REPEAT  LAST  STEP  BV  RESTORING  SAVED  INFORMATION. 

r  »**  +  *****  +  *  *.  *  *  ****  +  +  **  +  #***  +  ***  +  *>c***  +  ****  +  *  +  *  +  +  ******  +  *  +  **i('**  +  + 

160   IF  <NQ. EQ. NQOLD)  JSTART  -  1 


IF  <NQ. EQ. NQOLD)  JSTART  -  1 

T  -  TOLD 

NO  ■  NO OLD 

K  «  NQ  +  1 

GO  TO  1 20 
170  CALL  SUB1 <N, V, SAVE, CSAVE, VMAX, ERROR, PW, IP) 
170  CALL  LINK<'SUB1') 

230  CALL  SUB2CN, V, SAVE, CSAVE, VMAX, ERROR, PW, IP) 
230  CALL  LINK<  ' SUB2  '  ) 

GO  TO  <249, 798), JLNK 
240  QO  TO  <  250  ,  640  ),IRET 
250  CALL  SUB3(N, V, SAVE, CSAVE, VMAX, ERROR, PW, IP) 


+  H.  +  +  +  +  + 


250  CALL  LINKC  '  SUB3  '  ) 

**********  +  **  +  +  **  +  +  *.+  *  +  *i(.  +  *  +  +  *  +  +  +  +1>  +  *  +  *  +  +  *  +  +  *  +  **  +  +  +  +  +  *  +  +  +  ***  +  *  +  +  +  *'C  +  +  +* 

*  UP  TO  2  CORRECTOR  ITERATIONS  ARE  TAKEN.    CONVERGENCE  IS  TESTED  BV  + 

*  REQUIRING  THE  L2  NORM  OF  CHANGES  TO  BE  LESS  THAN  bND  WHICH  IS      + 
+    DEPENDENTONTHEERRORTESTCONSTANT.  * 

THE  SUM  OF  THE  CORRECTIONS  IS  ACCUMULATED  IN 

OR(I>.    IT  IS  EQUAL  TO  THE  K-TH  DERIVATIVE  0 

H**K'<FACTORIAL<K-1  )*A<K)  )  ,  AND  IS  THEREFOR 

Tur   orTimi  c o o  n  o  c   rn   Tur   inurer   nnur d   n  r   u   n 


ERRi 
BV 


IN  THE  ARRAV 

F  V  MULTIPLIED 
•rn\r>  i  /  ,  Hni'  i»  inLMLrunE  PROPORTIONAL 
THE  LOWEST  POWER  OF  H  PRESENT.   <H**K) 

•it**********************************-'*'-'' 


i*********- 


C 

c  + 

c 

c  + 

c* 

c* 

c* 

C*    TO  THE  ACTUAL  ERRORS  TO 

C  +  *  +  ***  +  4.*  +  4.  +  +  **  +  +  +  +  -4.  +  ****  +  ** 

DO  270  I  »  1  ,  N 
270     ERRORU)  -  0.0 

DO  430     L-1  ,  2 
C  CALL  DIFFUN<!T,  V,DV) 

CALL  L  I  N  K  <  '  DIFFUN') 

DO  280   1=1 , N 
280  C  S  A  V  E  <  I  ,  1  >  =>  DV<I) 

J  L  N  K  ■  1 

IF ( I WEVAL. LT. 1 )  JLNK  =  2 

IF'.  JLNK.EQ.2)  GO  TO  3  50 
C  +  +  ****  +  +  *  +  +  +  **  +  +  +  *  +  +  *  +  ****  +  ****  +  *  +  +  ****  +  *  +  +  +  +  **  +  +  +  **  +  *  +  +  +  ***  +  *  +  >i.  +  +  *  +  *  +  + 

C+  IF  THERE  HAS  BEEN  A  CHANGE  OF  ORDER  OR  THERE  HAS  BEEN  TROUBLE  * 

C*  WITH  CONVERGENCE,  PW  IS  RE-EVALUATED  PRIOR  TO  STARTING  THE  * 

C  *  CORRECTOR  ITERATION  IN  THE  CASE  OF  STIFF  METHODS.    IWEVAL  IS  * 

C*  THEN  SET  TO  -1  AS  AN  INDICATOR  THAT  IT  HAS  BEEN  DONE.  * 

310    DO  340  J  ■=  1  ,  N 
F  =  V< J,  1  ) 


330 
340 


R  ■=  EPS+AMAX1  <EPS, 
V<J,1)  =  V<J,1)  +  R 

D  =  A  <  1  )  *  H  y  R 
CALL  DIFFUN(T,V,DV) 
CALL  LINK< ' DIFFUN ' ) 
DO  3  3  0  I  =  1  ,  N 
PW ( I  +  <  J- 1  )*N  )  = 
V  <  J  ,  1  )  -  F 


ABS(F)  ) 


<DV<  I  )-CSAVE< 1,1) >*D 


35 

C  +  ****>!. +  *****#*#******#*♦*••*****••*****************•**  *********  41*41 4141414141 

C*    ADD    THE     H>ENTITV    MATRIX    TO    THE    JACOBIAN    AND    DECOMPOSE     INTO    LU    -    PH       * 

C#*  +  +  +  +**#****.*»  +  **  +  ***#«i**4i*»*«*i***  +  *****«i4l*****  +  fr*  +  *******  +  ****  + +  *  +  **  + 

296      DO  386  I  ■  1  H 
380      PW< I*(N+1 )-N>  -  1.0  ♦  PN<I*<N*1 >-M> 
IUEVRL  -  -1 
C  358  CALL  SUB7( N, V, SAVE, C8AVE , VMAX , ERROR , PH,  IP > 
330  COLL  L1NKC ' SUB7 ' > 

GO  TO  (448,410), JLNK 
C  410  CALL  SU84(N, V. SAVE, CSAVE, VMAX, ERROR, PH, IP) 
410  CALL  LINK( 'SUB4' > 

GO  TO  <430, 448), JLNK 
430  CONTINUE 
C  440  CPLL  SUB3(N, V , SAVE , Ct AVE , VNAX, ERROR, PH, IP) 
440  CPLL  LINKC 'SUBS' ) 

GO  TO  <?58, 715,558, 788, 549), JLNK 

C  +  +  * +  +  +  4*  +  *****  +  **  +  +  +  ****  +  +  **  +  ***#  +  «i4i*******-«**»*  +  #  +  *  +  +  +  *  +  +  ****  +  *  +  +  +  +  +  *  + 

C*  REDUCE  1HE  FAILURE  FLAG  COUNT  TO  CHECK  FOR  MULTIPLE  FAILURES.  * 

C*  RESTORE  T  TO  ITS  ORIGINAL  VRLUE  AND  TRV  AGAIN  UNLESS  THERE  HAVE  4. 

C*  THREE  FAILURES.   IN  THAT  CASE  THE  DERIVATIVES  ARE  ASSUMED  TO  HAVE  + 

C*  ACCUMULATED  ERRORS  80  A  RESTART  FROM  THE  CURRENT  VALUES  OF  V  IS  * 

C*  TRIED.  THIS  IS  CONTINUED  UNTIL  8UCCE68  OR  H  -  HMIN.  41 

Ct+***++*<+  +  *  +  *  +  +  +  +  *  +  +  »*»***  +  *****+»****  +  ***4.*  +  +  +  **  +  *  +  *  +  *  +  **+***  +  4.  +  *  +  ++* 

546   KFLAG  -  KFLAQ  -  2 

IF (ABS(H).LE.HMIN)  60  TO  748 
T  *     TOLD 

IF  (KFLAQ. LE. -5)  80  TO  720 
550  JLNK  -  2 
C  640  CALL  SUB6(N, V, SAVE, CSAVE.VMAX, ERROR, PH, IP) 
640  CALL  LINK( ' SUB4 ' ) 

GO  TO  (700,758, 170,250), JLNK 
700   DO  710  I  •  1  ,N 
710    VMAX(I)  -  AMAX1  (VMAXd),  A0S(V(  I,  1  >)> 

JSTART  -  NQ 
715  RETURN 
720   IF(NQ.GT. 1 )  80  TO  726 

IF(  ABS(H).LE.2.E8*HHIN>  80  TO  780 
GO  TO  550 
725  R  -  H'HOLD 

DO  738    I  -  1,N 
V( I , 1 >  -  SPVE(  1,1) 
738    V(I,2)  ■  SAVE<I,2)*R 
NQ  -  1 
KFLAG  >  1 
GO  TO  170 
740   KFLAG  ■  -1 
HNEH  -  H 
JSTART  -  NQ 
RETURN 
C  750  CALL  SUB8<N, V, SAVE, CSAVE, VMAX, ERROR, PH, IP) 
750  CALL  L  INK<  '  SUBS'  ) 

GO  TO  (  1  30,250,640,440), JLNK 
700   KFLAG  -  -4 
RETURN 
END 


Sample  user-provided  Overlay:  Y"  ■  -Y 


.TITLE  INTRO 

.GLOBL  N, RETURN, DIFFUN 

.CSECT  ADDR 

addr:  .*. 

.CSECT  COM 

T-. +74 

.  CSECT 

BEGIN:     MOV  ADDRH0.L1ST*! 
MOV  ADDR,LIST*4 
JSR  15, DIFFUN 
BR  RET 

LIST:      .WORD  T,0,8 

RET:       JSR  13, RETURN 

.END  BEGIN 


SUBROUTINE  D IFFUN( T , V , DV > 

DIMENSION  V(2,7>,  DV(2) 

DV( 1  >  -  V(2,  1  ) 

DV<2>  -  -V( 1,1) 

RETURN 

END 
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Regular  Overlays;  INTRO  vfth  thq  appropriate  aubroutine 
name  lotrodpcaa  eaqft  overly 


[NTHO 


.TITLE  IN.r 
.GL06L  N, RETURN 
.CSECT  ODDR 
ODDR:  .-. 
. CSECT 
BEGIN: 


,SUB1 


LIST! 

RET:      J 

.END  BEGIN 


I  15, SUB1 
BR  RET 

. word  n,  0, 0,0,0, 0, e,e 

JSR  15, RETURN 

M 


c  *  *  *  + 
C*  SE 
C*  TV 
C*  TH 
C*  BE 
C*  ST 
C*  SC 
C**** 


c  *  *  *  * 
C*  TH 


O 

c* 
c* 
c* 
c* 
c* 
c* 
c* 


OC 


sue 

COM 
1  I 
2ENQ 

DIM 

1 

***,* 
T  TH 
PE. 
IS  S 
CPUS 
EP  I 
OLIN 

*  *  *  * 
0(2 
IF 
GO 

*  *  *  * 
E  FO 
CUfift 


ROUTINE 
HON  'CO 
DOUB,  IN 

2.ENQ3, 
ENSION 
ER 

******* 

E  COEFF 

CHECK 
ECTION 
E  OF  TH 
F  IT  HP 
G  BEFOR 
******* 

>  *  -1  . 

(NQ  ,GT 
TO  (221 

******* 

LLONING 
CV  PERM 


SU61 (N, 

ns   r<7>, 

EVPL, IRE 
BND.BR.D 
V<N,7>,S 

ROR<N>  ,  P 

*  *  *  *  *  *  ** 


ICIE 
FOR 

SET 
E  OR 
S  NO 
E  EX 
*  *  *  * 
8 

.  6) 
,222 

COEFFIC 
ITTED  BY 


NTS 
EXCE 
IWE 
DER 
T  YE 
IT  I 
***  * 

NQ- 
,223 


V,BAVE, 
JLNK, JS 
T, IRET1 
EL, DELI 
0VE(N,7 
W(N,N), 
******* 

THAT  DE 
SSIVE  0 
YOL  .8T 
CHANGE, 
T  BEEN 
F  IT  HO 
******* 


C60VE, 
TORT.K 
,H,HOL 
,D 

>, esov 

I  P  <  N  > 

****** 

TERMIN 
RDER. 
.0  IF 

OND  T 
DONE  < 
S  BEEN 

****** 


YMOX, ERROR, PH, IP) 

,KFLOa,NQ,NQOLD,NEWQ, EPS, T, TOLD, 
D,HNEW,HMfiX,HMIN, E,EUP,EDUN, ENQ1 , 

E<N, 2) , YMOX(N), 


*****  ****************************** 
E  THE  ORDER  AND  THE  METHOD  * 

THE  LOST  TWO  STPTEMENTS  OF  + 
PW  IS  TO  BE  RE-EVOLUOTED  * 

HEN  REPEOT  THE  INTEGROTION  * 

IRET  •  1>  OR  SKIP  TO  0  FINAL       ♦ 

COMPLETED  <IRET  -  2).  * 

*********************************** 


,224, 223, 226), NQ 

******************** 

IENTS  SHOULD  BE  DEFI 
THE  MOCHINE.   THEV 


************************* 

NED  TO  THE  MOXIMUM 

ftRE,  IN  THE  ORDER  USED. . 


-1 
-2 

-  1 

-  1 
-1 

£  *  *  *  * 

221 
222 


223 


224 


225 


22* 


'3, -1 '3 
2/25,  -7/1 
20/274, -22 

90/44 1,-58 
********** 

ft  <  1  >  -  -  1 

GO  TO  230 

ft<  1  > 

P<3) 

GO  TO 

PC  1  > 

ft  ( 3 ) 

P(4> 

GO  TO 

0<  1  ) 

ft  (  3  ) 
.0(4) 

ft  < 5  > 

GO  TO 

ft<  1  > 


•1/5,-1/50 


*** 

* 
* 
* 
* 
* 
* 


3/274J-83/274, -15/274, -1x274 

/63,-15/36,-23/25  2,-3/232,-lx 

***************************** 


1  764 

*********+******+********. 


-0, 
-0, 

238 


238 


0(3) 
ft  <  4  ) 
0(5) 

ft<6> 
GO  TO 

ft  <:  i  >  - 

fi<3>  - 

ft  <  4  >  - 

0(5)  = 

0(6)  » 

ft<7>  • 

RETURN 
END 


ft( 
-8 

230 
-0 
-0 

-0 
-0 

230 
-0 
-0 
-8 
-8 
-0 

230 
-8 

-0 

-8 

-8 

-0 
-8 


809808888 

6666666666666667 
3333333333333333 

8.  5454543454545453 

1  > 

. 89898909098909891 

. 480088088 
. 708008800 

.280  8  00000 
. 020008080 

. 437956204379562 
. 821 1 6790321 1 6788 
.3  102189781821898 
. 05474452S54744526 
. 8836496350364963384 

.4861632653861 225 
. 9206349206349286 
.  4 1  * 666 66 6 666 66 6 7 
. 0992863492863492 
.81 1 984761 984761 9 
.008566893424036282 
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SUBROUTINE  8UB2 < N , V , 8AVE , C8 AVE , VMAX , ERROR , PW ,  IP) 

COMMON  •COM/'  A<?>,  JUNK,  JST ART , K. KFL AQ ,  NQ , NQOLD , NEWQ ,  EPS,  T,  TOLD, 
1   IDOUB, IWEVAL, I  RET, IRET1 , H, HOLD, HNEW, HMAX, HM I N, E , EUP, EDWN, ENQ1 , 
2ENQ2,ENQ3,BND,BR,DEL,DEL1 ,  D 

DIMENSION  V<N,7),SAVE<N,7>,CSAVE<N,2>, VMAXCN), 
1  ERROR<N>,PM<Nf N>, IP<N> 

DIMENSION  PERTST<7,3> 

C*  +  +  +  +  *  +  +  +  *  +  +  +  +  +  *  +  +  +  +  +  **m  +  **  +  *<*  +  +  +  +-***  +  +  +  +  +  +  m  +  +  +  *  +  +  *  +  **.  +  ±  +  +  +  **  +  +  4,*  +  +  +  t,  +  4, 

C*  THE  COEFFICIENTS  IN  PERT8T  ARE  USED  IN  8ELECTINQ  THE  STEP  AND         * 
C*  ORDER,  THEREFORE  ONLV  ABOUT  ONE  PERCENT  PCCURRCV  IS  NEEDED.  * 

DATA  PERTST/2.,4.5,7.333,10.42,13,7,17.15,1., 

1  3. ,6. , 9. 1 67, 12.5, 15.98, 1 . , 1 . , 

2  I  .,  1  .,.  5,  .  1  667,  .  94U7,  .  808333-' 
JLNK  »  1 

230   K  -  NCi+l 
IDOUB  -  K 

ENQ2  ■  . 5/FL0AT<NQ  +  1 > 
ENG3  -  .5/FL0AT<NQ  «■  2> 
ENQ1  -  0.5/FLOAT<N8> 
EUP  -  <PERTST<NQ,2>*CPS>**2 
E  -  <PERT6T(NQ, 1 >*EPS>**2 
EDWN  -  <PERTST<NQ, 3)*EPS)**2 
IF(EDHN. EQ. 0)  JLNK  -  2 
OND  -  <EPS*ENQ3)**2 
240  IWEVQL  -  2 
RETURN 
END 


SUBROUTINE  SUB3<  N  ,  V  ,  SAVE  ,  C8AVE  ,  VMAX  ,  ERROR  ,  PW  ,  I  P  "> 

COMMON  'COM'  A<7) , JLNK, JSTRRT, K, KFLAQ, NQ, NQOLD, NEWQ, EPS, T, TOLD, 

1    IDOUB, IWEVAL, IRET.IRET1  , H , HOLD , HNEW , HMAX , HM I N , E , EUP , EDWN , ENQ 1  , 

2ENQ2,ENQ3,BND,BR,DEL,DEL1  ,  D 
DIMENSION  V(N,7),SAVE<N,7),C8AVE<N,2),VMAX<N), 

1  ERRORCN)  , PW<N,N> ,  IP<N> 

C******+********. ***+******•*******************************♦.*** **♦+**••« 

C*  THIS  SECTION  COMPUTES  THE  PREDICTED  VALUES  BV  EFFECTIVELY  * 

C*  MULTIPLVINQ  THE  SAVED  INFORMATION  BV  THE  PASCAL  TRIANGLE  * 

C*  MATRIX.  * 

£••*******+*********•**+******•*******+*****+*******+*******+♦+*+*++*•** 

250   T  -  T  ♦  H 

DO  260  J  -  2,K 
DO  268  J1  -  J,K 

J2  ■  K  -  J1  ♦  J    -     1 
DO  260  I  -  1  ,  N 
260         V<I,J2)  -  V<I,J2)  ♦  V<I,J2+1) 
RETURN 
END 


SUBROU1 

COMMON 

1    IDOUB 

2EN02.EN 

DIMENSI 

1 

C************ 

C*  CORRECT  AN 

C*  ESTIMATE  T 

C*  BR*2*DEL  A 

C*  IS  SAID  TO 


INE  SUB4<N, V, SAVE, C6AVE, VMAX, ERROR, PW,  IP> 

/COM/  A< 7 >  JLNK, J8TART,K,KFLQG, NO, NQOLD, NEWQ, EPS, T, TOLD, 

,  IWEVAL,  I  RET,  1  RET  1  , H , HOLD , HNEW , HMAX , HM  I  N , E , EUP , EDWN , ENQ 1  , 

03,  0ND.BR,  DEL,  DELI  ,  0  «u,M<k 

ON  V<N,7),SAVE<N,7),CSAVE<N,2>,VMAX<N), 

•S2222i2iiS2i2i2iiISi2i 

D  COMPARE  DEL,  THE  L2  NORM  OF  CHANGE /VMAX ,  WITH  BND.        * 

Se  5SlSe  OF  THE  L2  NORM  OF  THE  NEXT  CORRECTION  BV  • 

ND  COMPARE  WITH  BND.    IF  EITHER  IS  LE*S,  THE  CORRECTOR      * 

HAVE  CONVERGED.  ^- ** *****••• 


418 


420 

426 


DEL  -  0 
DO  42 

V<  I 
V<  I 
ERR 
DLL  -  D 
CON 
IF<L.GE 
DELI  - 
JLNK  - 
IF(AMIN 
RETURN 
END 


.  0E0 
0  I 
,  1  > 
.2) 
0R<  I 
EL 

TINUI 
.  2> 
DEL 

1 <DEL,BR*DCL*2.0>   .LE.   BND)   JLNK  -  2 


-  V<  I,  1  )  ♦  A<1 >*CSAVE<I,  1  > 

-  V< I , 2)  -  CSAVE<  1,1) 

)  -  ERR0R<  I )  ♦  CSAVE<  1,1) 
<  C  S  A  V  E  <  I  ,  1  >/VMAXU)>**2 

E 


JE 
BR 


-  ANAX1 < . 9*BR,  DEL/DEL1) 
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c  *  *  +  * 

C+  TH 
C+  PO 
O.  TH 
C*  MA 

C  *  IS 
C*  ST 
C  *  *  +  + 

440 


468 
470 
480 


SUBROUTI 

COMMON  / 

t    IDOUB, 

2ENQ2, ENQ 
DIMENS  10 
1 

********* 

E  CORRECT 
SSIBILITI 
IS  IS  EIT 
TRIX  PW  H 
TAKEN.  0 
EP  IS  RED 
********* 

GO  TO  (4 
T  =  TOLD 
I T  <  QBS  C  H 
IF< IHEVP 
IWEVPL  - 
I  RET  1  - 
JLNK  -  1 
GO  TO  33 
KFLPG  ■ 
NO  -  NGO 
DO  480  I 
DO  460 
V<I, 
H  -  HOLD 
JSTART  » 
JLNK  •  2 
GO  TO  53 


NE  6 
COM' 
IWEV 
3,BN 
N  V< 
ERRO 

*  *** 
OR  I 
ES  P 
HER 
AS  P 
THER 
UCED 

*  +  *  * 
40,  4 

>.LE 

L.  NE 

2 
2 

6 

-3 

Lt> 
-  1 
J  - 

J)  - 

NO 


UB5<N .V.eflVE.CSPVE.VMPW. ERROR, PW, IP) 

A<7), JLNK, J8T ART, K.KFLA8, NO. NQOLD.NEWQ, EPS, T, TOLD, 
AL.IRET.IRETI.H.HOLD.HNEW.HMAX.MMIN.E.EUP.EDWN.ENQI, 
D, BR, DEL, DELI  , D 

N,7),SflVE<N,7>,C8flVECN,2),VMPK<N), 
R<N>, PN<N,N>, IP<N> 

************»**********•***********+*»**************** 
TERPTJON  FAILED  TO  CONVERSE  IN  2  TRIES.  VARIOUS  * 
RE  CHECKED  FOR.    IF  H  19  ALREADV  HMIN  AND  * 

ADAMS  METHOD  OR  THE  STIFF  METHOD  IN  WHICH  THE  * 

LREADV  BEEN  RE  -EVALUATED ,  A  NO  CONVERGENCE  EXIT  * 
WISE  THE  MATRIX  PW  IS  RE-EVALUATED  AND'OR  THE  * 

TO  TRV  AND  GET  CONVERGENCE.  * 

a***************************************************** 

90,470, 470), JLNK 


.HMIN  .AND.  <  IWEVAL-1  )  .LT. -1  >  GO  TO  463 
.0)  H  -  H*0,2SE0 


,N 
1  ,K 

SOVC<I, J> 


C*  THE  CORRECTOR  CONVERSED  AND  CONTROL  16  PASSED  TO  9TATEMENT  320  • 

C*  IF  THE  ERROR  TE6T  19  O.K.,   AND  TO  540  OTHERWISE.  * 

C*  IF  THE  STEP  IS  O.K.   IT  19  ACCEPTED.  IF  IDOUB  HAS  BEEN  REDUCED  * 

C*  TO  ONE,    A  TEST  IS  MADE  TO  8EE  IF  TME  9TEP  CAN  BE  INCREASED  * 

C*  PT  THE  CURRENT  ORDER  OR  BV  GOING  TO  ONE  HIGHER  OR  ONE  LOWER.  • 

C*  SUCH  H  CHAN8E  IS  ONLV  MADE  IF  THE  STEP  CAN  BE  INCREASED  BV  AT  * 

C*  LEAST  1.1.    IF  NO  CHANGE  IB  POSSIBLE  IDOUB  IS  SET  TO  8  TO  ♦ 

C*  PREVENT  FUTHER  TESTING  FOR  9  STEPS.  * 

C*  IF  P  CHAN6E  IS  POSSIBLE,   IT  IS  MADE  AND  IDOUB  IS  SET  TO  • 

C*  NQ  ♦  1      TO  PREVENT  FURTHER  TESING  FOR  THAT  NUMBER  OF  STEPS.  • 

C*  IF  THE  ERROR  HAS  TOO  LARGE,  THE  OPTIMUM  STEP  SIZE  FOR  THIS  OR  * 

C*  LOHER  ORDER  IS  COMPUTED,  AND  THE  9TEP  RETRIED.    IF  IT  SHOULD  * 

C*  FAIL  TWICE  MORE  IT  IS  AN  INDICATION  THAT  THE  DERIVOTIVES  THAT  * 

C*  HAVE  ACCUMULATED  IN  THE  V  ARRAV  HAVE  ERRORS  OF  THE  WRONG  ORDER  * 

C*  SO  THE  FIRST  DERIVATIVES  ARE  RECOMPUTED  AND  THE  ORDER  IS  SET  * 

C*  TO  1  .  * 


1  ,N 
(.  ERROR*  I  >/VMAX<  I  >  >**2 


496   D  »  0.0 
DO  580  I 

508      D  ■  D 

IWEVPL  -  0 
IF(D.GT.E)  JLNK  *  5 
IF< JLNK. EQ. 5)  GO  TO  536 
IF  <K. LT . 3)  GO  TO  520 


C*  COMPLETE  THE  CORRECTION  OF  THE  HIGHER  ORDER  DERIVATIVES  AFTER  P 
C*  SUCCESFUL  STEP. 


51  0 
520 


5  30 

5  35 
536 


3,K 
-  1 


N 
V<I, 


DO  510  J  - 
DO  510  I 
V( I , J;  = 
K  F  I.  P  G  .=  *  1 
H  N  E  W  »  H 
IFCIDOUB.LE.  1 
IFCJLNK.EQ.3) 
IDOUB  -  IDOUB  ■   1 
IF(IDOUB.GT.I)  GO  TO 
DO  530  I  =  1  ,  N 

CSPVE^ 1,2)  -  ERRORC  I  > 
JLNK  =  4 
RETURN 
END 


j>  ♦  r<  j:>*error<  i  > 


)  JLNK-3 
GO  TO  536 


535 
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SUBROUTINE  6U86<N  ,  V  ,  SAVE  ,  CSflVE  ,  VMfiX  ,  ERROR  ,  PW  ,  IP> 

COMMON  'COM'  fl<7>.JLNK.  J«T«RT , K , KFLhO, NQ , NQOLD , NEWQ ,  EPS,  T,  TOLD, 

1    IDOUB,  IWEVfiL,  I  RET, IRET1 , H, HOLD , HNEW , HMfiX , HM I N , E,EUP,EDWN,ENQ1, 

2ENQ2,ENU3,BND,BR,DEL,  DELI  ,  D 


,  *■  u  n  u  u  l  v     dc  u  i  v  i  u  c  v     hi   vn 

C+  ONE  HIGHER  RESPECT  I VELV . 

GO  TO  (680, 358, 640) , JLNK 
550   PR2  -  <D'E)»*ENQ2*1  .8 
PR3  «  1 .E*28 

IF<<NQ. GE.6).0R. <KFLfi8.LE.-1 >>80  TO  576 
r\     ^     a     a 


560 


0  =  0.0 

DO  560  I  -  1  ,N 

0  ♦  <<ERROR<I)  -  CSAVE< I, 2>>'VMfiX<I>>**2 


PR3  «  <D'EUP)**ENQ3*1 . 4 

PR1  =  1  . E  +  20 

IF  <NQ. LE.  1  )  00  TO  590 

D  -  0.  0 

DO  58  0  I  -  1  ,  N 

0    *    D  ♦  <V<  I  ,  K)/VMQX<  I  >  >«"»2 
PR1  -  <D/'EDWN)**ENQ1  *1  .  3 
CONTINUE 

IF  <PR2,LE.PR3>  00  TO  650 
IF  <PR3. LT. PR1 >  80  TO  660 
R  -  1 .©^flMfiXI <PR1 , 1 . E-4> 
NEWQ  *  NQ  -  1 
610  IDOUB  -  8 

IF  ( <KFLfiG. EO. 1 ) .OND. <R. LT. < 1 . 1 > ) >  60  TO  694 
IK  <NEWQ.LE.NQ)  00  TO  630 


570 


580 
598 


680 


C*  COMPUTE  ONE  ADDITIONAL  SCALED  DERIVATIVE  IF  ORDER  IS  INCREASED 

i+  +  ***  +  +  +  +  +  +  *  +  *  +  +  +  +  +  +  +  +  **  +  +  *+  +  +  +  *  +  +  **+  +  *  +  +  +  +  *+  +  *  +  *   +   *  +  *  +  **  +  +<+  +  +  *H 

DO  620  I  -  1  ,  N 


628 
630 


640 


658 


660 


670 


680 


690 

6*4 
695 


DO  620  I  -  1  ,  N 

V<I,NEWG*1>  -  ERROR* I >*A<K) 'FLOAT(K) 
K  »  NEWQ  +  1 

IF  (KFLAG.EQ. 1 )  80  TO  670 
H  -  H*R 
IRET1  -  3 
JLNK  •  2 
GO  TO  695 

IF  (NEWQ. EO. NO)  JLNK-4 
IF<JLNK.EQ.4>  GO  TO  695 
NO  ■  NEMQ 
JLNK  »  3 
GO  TO  695 

IF  <PR2. QT. PR1 >  00  TO  600 
NEWQ  -  NO 

R  -  1 . 0'HMfiX1 <PR2, 1 . E-4> 
GO  TO  6  10 

R  -  1.8'HMAX1tPR3,1.E-4) 
NEWQ  -  NQ  ♦  1 
'3  0  TO  618 
IRET  -  2 

R  ■  AMIN1  <R, HMAX/ABSCW)  > 
H  -  H*R 
KNEW  -  H 

IF  «:Nw.  EO.NEWO)  GO  TO  688 
NQ  -  NEWQ 
JLNK  -  3 
GO  TO  695 
R1   »  1.0 
DO  *  9 8  J  -  2 ,  K 

R1   -  R1 *R 

DO  6  9  8  I  -  1 , N 

V  <  I  ,  J  >  -  V  <  I  ,  J  >  *  R  1 
IDOUB  =  K 
JLNK  -  1 

RETURN 
END 


++**+++ 


OUTPUT. 

fl<  I  ,  J 

P  <  I  ,  J 

IP(K> 

IP<N> 

USE  'SO 

DETERM< 

IF  IPCN 


iRDER  OF  MPTRIX 

DECLARED  DIMENSION  OF  PRftftV  P. 
IPTRIX  TO  BE  TRIANGULRRJZED.  <FOR  8TIFF  METHODS, 
precision;   ALL  OTHER  VfiRIQBLE  PRE  DOUBLE 


p  is  single 
precision. > 


I 


LE.J   »  UPPER  TRIANGULAR  FACTOR,  U. 

I.QT.J    -  MULTIPLIERS   •  LOWER  TRIPNQULPR  FACTOR,  I-L, 
K.LT.N  »  INPEX  OF  K-TH  PIVOT  ROW, 
<-1 >**<NUMBER  OF  INTERCHANGES)  OR  0. 
LVE'  TO  OBTAIN  SOLUTION  OF  LINEPR  SVSTEM. 
P>  -  IP<N>*fl<1,1>*fl<2,2)*.,,*fl<N,N>, 
>  «=  0,  A  IS  SINGULAR,  SOLVE  WILL  DIVIDE  BV  ZERO, 


f  01 


TO  5 


■IP<N> 


350 
360 


DIMENSION  A<N, N> ,B< 1 > , IP<N>, V< 1 > 

GO  TO  ■;  1  0  1  ,  3  5  0  )  ,  J  L  N  K 

J  L  N  K  =  2 

N  D  I  M  =  N 

I  P  <  N  >  -  1 

DO  6   K=1  ,  N 

I  F  <  K  .  E  6  .  N  >  GO 

KP1  -  K+1 

U   m       U 

DO  1    I-KP1,N 

IF<RBS<P<I,K)>.QT.PBS<P<M,K>>> 
CONTINUE 
IP<K>  -  M 
IF(M.NE.K)    IP<N> 
T  -  fl<M, K> 
P  <  M  ,  K  )  -  P  t  K  ,  K  > 
fl  <  K  ,  K  >  •»     T 
I  F  <  T  .  E  Q  .  0  >  GO  TO 
DO  2   I = K  P 1 
P< I , K>  ■ 
DO  4   J-KP1 , N 
T  -  fl  <  M  ,  J  > 
fl<M, J)  =  P<K 
fl  <  K  ,  J  )  ■=  T 
IFCT. EQ. 0.  ) 
DO  3     I-KP1 , N 

P< I , J)  =  P<  I  , J>  ♦  R(| ,K>*T 
CONTINUE 
IF<A<K,K>.EQ.0 
CONTINUE 

IF< IP<N> . EO. 0>  JLNK  ■  1 
IF< JLNK. EQ.  1  > 
360  Ixl.N 


M-l 


N 


-P< I , K>  /  T 

'<i 

,  J> 

GO  TO  4 

P< I, J> 

IP(N> 

IF<IP(N).EO.0> 

IFC JLNK. EO. 1 )  GO  TO  10 

DO 


0 


B<  I  > 

JLNK 


V< I+N>  -  BC  I  >*H 


SOLUTION  OF  LINEAR  SVSTEM,   A*X  -  B. 

INPUT. . . 

N  -  ORDER  OF  MPTRIX. 

NDIM  -  DECLARED  DIMENSION  OF  PRRPV  P. 
fl  =•  TRIANGULARIZED  MATRIX  OBTAINED  FROM 
6  -  RIGHT  HAND  SIDE  VECTOR. 
IP  =  PIVOT  VECTOR  OBTAINED  FROM  'DECOMP 
OUTPUT. . . 

B  -  SOLUTION  VECTOR,  X. 
I F  <  N . E 6 .  1  >  GO  TO  « 
NM1  =  N-1 
DO  7   K -  1  , N  M 1 

KP1  =  K  +  1 

M  -  IP<K> 
T  =  B  <  M  > 


K)*T 


K  =  K  M 1   +  1 

B  <  K  >  =  B  <  K  >  /  fl  (  K  ,  K  > 

T  =  -  6  <  K  > 

00  8   1=1 ,KM1 
*  B  <  I  >  -  B  <  I  >  ♦  A  <  I  ,  K 
">  B  <  1  >  =  B  <  1  )  /  P  <  1  .  1  > 
10  RETURN 

END 


DECOMP 


B(M  ■■> 

-    B<K> 

B  <  K  ) 

=     T 

DO     7 

I-KP1 , 

N 

B<  I  •> 

=    BCD 

+ 

P<  I 

DO   '8 

KB-1  , 

NM1 

KM1     = 

:     N     -     KB 

i 

*T 


ill 


C*  TH 
C  *  TO 

C  +  *•  *  * 
759 


76S 
770 

76* 
783 


SUBROUTINE  6UB6 < N , V , SAVE , CS AVE , VM AH , ERROR , PW , IP > 

COMMON  -'COM/  A  <  7  >  ,  JLNK  ,  J8TART  ,  K  ,  KFLfiS  ,  NO.  ,  NQOLD  ,  NEWO  ,  EPS  ,  T  ,  TOLD  , 

1    IC)OUB,IWEVPLIIRET>!RETl,M,H01.D,HNEW,HMPX)HMIN,EIEUP)Et'WN,ENi?1, 
2ENQ2,ENQ3,BND,BR,  DEL,  DEL  1  ,  t> 

DIMENSION  V(N,7),S«VE<N,7),CS0VE<N,2),VM0X(N), 
1  ERROR'N) ,PW<N,N>, IP<N> 

*V++++>**+4V  +  +  +  *  *  +  ***  +  *  +  **-K*  +  ******  +  ***  +  **t.+  +  +  +  +  ****  +  4'*  +  +  +  +  +  +  +  +  +-l<*  + 

IS  SECTION  SCALES  ALL  VARIABLES  CONNECTED  WITH  H  AND  RETURNS        * 
THE  ENTERING  SECTION.  * 

H  =  RMfiXI  <HMIN, RMIN1  <  H  , HMAX)  > 
R1  =  1.0 
DO  7<i0  J  -  2,K 
R  -  H'HOLD 
R1  -  R1 *R 

DO  760  I  *  1  ,N 

V<I,J>  -  SAVE< I  , J>*R1 
DO  7  7  0  I  -  1  , N 

V  <  I  ,  1)  »  S  R  V  E  <  I  ,  1  > 
IDOUB  -  K 
JLNK  »  I  R  E  T  1 
GO  TO  7  8  5 
K.KLAG  -  -4 
JLNK  -  4 

RETURN 
END 
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